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Guido Brunnett 
Department of Mathematics 
Naval Postgraduate School 

Peter E. Crouch 

Center for Systems Science and Engineering 
Arizona State University 

January 28, 1993 



Abstract 

This paper deals with the derivation of equations suitable for the 
computation of elastic curves on the sphere. To this end equations 
for the main invariants of spherical elastic curves are given. A new 
method for solving geometrically constraint differential equations is 
used to compute the curves for given initial values. A classification of 
the fundamental forms of the curves is presented. 



1 Introduction 

Two current trends in Geometric Modeling are concerned with 

• the development of spline techniques on surfaces (see [14], [12], [10]) 

• the use of nonlinear spline curves of minimal elastic energy for the 
modeling of smooth shapes (see [2], [3], [4]). 

This paper blends both topics in considering elastic curves on the sphere. 

In euclidean space an clastic curve can be viewed as an arc length para- 
metrized ’’cubic” spline in tension, i.e. an elastic curve is a critical point of 
the functional ^ 

^(y) = / < y",y" > +(T < y\y' > ds .(1) 

Jo 
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in the space F of smooth maps 

j,:(0,/]-R”, |!/'l = l 

y{0) = Po, y(l) = Pu !/'(0) = V„, y'(l) = V, 

where Pq, Pi G R”, Vq G Tp^BP, Vi G Tp^BP, cr ^ R are fixed and / is 
variable. Expressing y” in the Frenet frame yields the functional $ in the 
form 

$(y)= f{K{s)f + ads (2) 

Jo 

which explains the interest in elastic curves as those curves that minimize 
bending. 

The notion of cubic splines can be generalized to curves on a Riemannian 
manifold M by replacing the usual derivative of the tangent vector field 
y' by the covariant derivative compatible with the metric of M (see [13]). 
Generalizing the functional (1) in this way, one obtains the concept of elastic 
curves on arbitrary manifolds. In the case of surfaces embedded in R^ the 
algebraic value of the covariant derivative of the tangent vector field ij' of a 
surface curve y is called the geodesic curvature Kg of y (see [5] and section 
2). Therefore we may define an elastic curve on a surface S' : /I C R^ R^ 
as an extremal point of the functional 

^(2/) = / (s(^))' (3) 

Jo 

in the subset F of F formed by curves lying on S. 

In section 3 a set of differential equations for elastic curves on the sphere is 
derived. This set includes a differential equation for the geodesic curvature of 
spherical elastica. Since the normal curvature of a spherical curve is constant, 
the differential equation for the geodesic curvature suffices to compute the 
ordinary curvature of a spherical elastica. Furthermore, a formula is given 
that expresses the squared torsion of a spherical elastica as a rational function 
of its curvature. 

In section 4 we describe the numerical algorithm used to integrate the 
set of differential equations derived in section 3. The equations have a very 
particular structure defined by a number of constants of motion, and in 
particular they constrain the elastic curves to lie on the sphere. We employ 
an algorithm introduced by Crouch and Grossman [8] which preserves the 
constraints exactly. 

Since Euler’s fundamental work on plane elastic curves it is known that 
these curves can be classified according to their shape. The tools developed 
in this paper enable us to present the fundamental forms of spherical elastica 
in the last section. 
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2 Geometric Preliminaries 



Let S : U C R^ denote a regular parametric surface and N the unit 

normal vector field of S. A curve a; ; / C R — + R^ is a curve on the surface 
S if and only if a; = 5 o c where c :/—♦(/ is a plane curve in U. The unit 



normal of S along a surface curve x will be denoted by n ;= o c. 

The Darboux frame 61,62,63 along x is the orthogonal frame defined by 



with u{t) = |a:'(i)|. 

The functions Kg. and Tg are called geodesic curvature, normal cur- 
vature and geodesic torsion. The geodesic torsion and the absolute value of 
geodesic and normal curvature are invariant under reparametrization of the 
surface. 

The geodesic curvature of a surface curve a; at a point x{i) is the ordinary 
curvature of the plane curve generated by orthogonal projection of x onto 
the tangent plane of S at x(t). It can be computed using the formula 



A surface curve with identically vanishing geodesic curvature is called a 
geodesic of the surface. 

The absolute value of the normal curvature of z at a point x{t) is the 
curvature of the intersection of S with the plane through x{t) spanned by 
the vectors x'{t) and n{t). While the geodesic curvature is the curvature of a 
surface curve from a viewpoint in the surface, normal curvature measures the 
curvature of the curve that is due to the curvature of the underlying surface. 
If K denotes the ordinary curvature of the space curve z the identity 




The equations that express the derivatives 6'i, 62, 63 in the Darboux basis 
61, 62, 63 are given by: 




(4) 

(5) 

( 6 ) 



K^ = K^ + K 



.2 

'3 



( 8 ) 



holds. 
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The geodesic torsion of a surface curve x at a point x{t) is the torsion 
of the geodesic that meets x at x(i) with common tangent direction. A 
curvature line of x, i.e. a curve with a tangent vector that points into one of 
the principal directions of the surface, is characterized by vanishing geodesic 
torsion. 



3 The differential equations of spherical elas- 
tica 



Let 5 be a parametrization of a patch on a sphere S 2 C of radius r and 
center 0 such that 

S = rN 

and X be an arc length parametrized (i.e. |x'| = 1) curve on S. Then 

f t V 

X = rn = ro^. 



In this situation (6) implies that 

Kn — , Tg — 0 . 

r 



(9) 



Since the absolute value of and Tg are invariant under reparametrization, 
these equations imply that any spherical curve is a curvature line with con- 
stant normal curvature. 

We consider the variational problem of minimizing 




“h (7 ds 



in the space F of C°° smooth maps 



y : |0,/1 ^ Sj, |i/'| = 1 

y{0) = Po, !/(l) = P,. /(O) = Vo, y'(l) = V, 

where Pq, Pi € S 2 , Vq € Tp^So, Vi 6 Tp^S 2 , cr € R are fixed and 1 is 
variable. 

From (4) one obtains the relation 



so that we wish to minimize the functional 



where 




under the constraints 



|ii|2 = l, b^ = x\ |x|^ = r^ 



Hence, we can apply the Euler- Lagrange equations to the functional 
F = |6;p + 6 + A(|6i| 2 - 1) + ^(|xp - r2) + 2 < A,x' - 6i > 
to obtain the differential equations which govern the extremals: 





0 

II 

1 

H 


(10) 




< 

II 

1 


(11) 


Combining these equations yields 






X'b, + xb[ - b'C = fix. 


(12) 


The derivatives of bx 


expressed in the Darboux basis are given by: 




K 


= Kgb2 bj 

r 


(13) 


K 


= - («J + ^)^l 


(14) 


K 


= {—3KgKg)bi + {Kg — Kg — 


(15) 



Substituting these derivatives and x = rb^ into (12) and rearranging gives 



(A' -f 3KgKg)bl 

+ “ '^9 + ^5 + 

- (7 + /^^ + ;(«? + ^))^3 = o. 

Finally, the linear independence of the vectors 61 , 62 ^.nd 63 implies that 

A = “ 2^5 ^ 

and 

Kg + -K^g — (C -f = 0. (17) 
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In order to determine the constant C in terms of the tension parameter 
cr, we consider the boundary condition 

for the extremal x. This condition is implied by the fact that the total length 
I of the curve is variable in the variation (see, e.g., [1]). Thus, 

-'=?(') -p + «-2<A(0,i'(')>=0. (IS) 

Substituting A according to (11) into the scalar product < A, x' > yields 
<A(/),x'(/)> = <A 6 i(/), 6 i(/)>-< 6 '/(/), 6 a( 0 > 

= - 5 («,('))= + ;^ + C. 

Substituting this expression for the scalar product into ( 18 ) we obtain 

C + = (19) 

These results are summarized in the following theorem. 

Theorem 1 An elastic curve x under tension a on the sphere of radius r 
satifies the differential equations 



x[ \ 




/ 


0 


1 /r 


0 \ 


/ Xl ^ 


X2 


= 




— l/r 


0 




X 2 


^3 / 




1 


0 


~^g 


0 y 


V ^3 y 



where X\ = x, X2 = rx' ^ X3 = x X x' and where the geodesic curvature 
X is a solution of 



( 20 ) 
Kg of 
( 21 ) 



The curvature of a spherical elastic curve can be obtained from (8) and the 
fact that the normal curvature of spherical curves is constant. The squared 
torsion of a spherical elcistica can be expressed as a rational function of its 
curvature. 



Theorem 2 The curvature k and the torsion t of a spherical elastic curve 
obey the relation: 



r T 



1 . 1 



1 .3 



k'* = + -?(t - o ) + ^1- 



w r 
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Proof: For the invariants k and r of an arc length parametrized curve x in 
the relation 

rK^ = |i.,,6;,6"l (22) 

holds where [a,b,c] denotes the determinant of three vectors in 

For a spherical curve b[ and b” can be expressed in the Darboux basis as 
follows: 



b[ — Hgb^ 63 

r 

b'l = -(— + K])bi + K'gb-2. 

Substituting the derivatives into (22) and squaring yields 






Since the differential equation (21) can be integrated to 



1 



1 cr. 



one obtains the claimed equation using (8). 



4 Tracking elastic curves on the sphere 

The problem we consider here is that of numerically integrating equations 
(20) and (21) in Theorem 1. One can of course simply integrate the equa- 
tions using a standard numerical package, such as an IMSL Runge Kutta 
routine. However, the system of equations possesses a very special structure. 
As pointed out in [3], equation (21) may be Integrated direct!}' in terms of 
Jacobi’s elliptic functions. We give more details of this process in the next 
section where we classify the various extremals. As for equation (20) we 
note that the components of the state vector [a:f , 1^, I3 satisfy algebraic 
constraints consistent with the fact that the matrix [a:i,X2,a:3] is simply a 
multiple r, of a rotation matrix. When a standard integration package is 
applied to the set of differential equations (20) and (21), these constraints 
are not preserved exactly, and in particular the norm of the vector Xi will 
not remain at the constant value This is a particularly important fact 
when we wish to integrate the equations over a large number of time steps 
and visualize the resulting curves. 

We have therefore made use of a new class of integration algorithms de- 
veloped by Crouch and Grossman [8], and Crouch, Yan and Grossman [7] 
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wliich do indeed preserve such structures. The algorithms are therefore called 
geometrically exact [6]. We briefly indicate the important aspects of these 
geometrically stable algorithms which pertain to the equations (20) and (21). 
Suppose that we wish to numerically integrate an ordinary differential equa- 
tion on R” given by the equations: 

x{t) = F{i,x{l)), 3- e R”, a:(0) = xo (23) 

where 

N 

F{t,x) = J2a^{t,x)Aj{x) (24) 

j=i 

n > Aj are vector fields on and are functions on R X R”'. Suppose 
that in addition we are given a set of functions on R’^ whose numerical values 
are constant along solutions of the equation (23) and that the level sets of 
these functions are manifolds. Denote the level set through Xo by AL It is 
convenient to assume the slightly stronger assumption that the vector fields 
Aj are everywhere tangent to A7. We also assume that there is an oracle that 
can integrate any vector field of the following form, to any desired accuracy: 

Z(x) = f;aM,(x). (25) 

j=l 

Here are real numbers. We define vector fields by setting: 

J=1 

and note that is simply the vector field F with coefficients '"frozen” at 
the point p. If we denote the flow of any vector field Z by (/,x) — > 6z{t^x)^ 
R X R^^ — > R””, then since the vector fields Aj are everywhere tangent to A/, 
it follows that x ^ M implies that 0pp{i^x) G M for all p and for all t for 
which the flow is defined. 

We now introduce the (explicit) geometrically exact Runge Kutta algo- 
rithms as described in [8]. Let Xk = x{tk) be a point of the integral curve 
X of (23). Then, define vector fields on R”” by freezing coefficients of F at 
various points as follows: 



N 



Cl (a:) = ^a^(4,a:fc)/ij(x) 
j=i 



N 

+ hc2l,0Fi{hc2l,Xk))Aj{x) 



J=l 
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^ 3 ( 3 ;) = + h{c3i + C32),0f^{hc32,0p^{hc3i,p)))Aj{x), 

i=i 

etc. Second, we define the numerical integration algorithm via an update 
rule: 

^/:+I — ^Fr (^^r! ^Fr_i (^^r— 1 i ■ ■ ■ ! ^Fi (^^1 1 ))) (26) 

where h is the “step length” and c,- and c,j are constants to be determined. 
'These constants are determined from the “consistency equations,” obtained 
by making the Taylor expansions in h, about h = 0, of both sides of (26), 
using on the left hand side the expression = 0p{h,Xk)- If the coefficients 
of /d agree up to i = </, then q is said to be Older of the resulting algorithm. 
Note that in general we have r > q, while for classical Runge Kutta schemes 
we' can always take r = q. Note that the update rule defined by equation 
(26) has the property that if a’*,, lies in M , then so does since each flow is 
defined by a vector field F with frozen coefficients. In the special case where 
n = jV, A1 = R’’‘ and A, = is the standard fth basis vector in /?”, the 
algorithm reduces to the form of a classical explicit Runge-Kutta algorithm. 

In the paper [8] the consistenc}^ equations are derived for the geometri- 
cally exact Runge Kutta algorithms via a careful geometric analysis of the 
equations (26). The results show that a geometrically exact third order ex- 
plicit Runge Kutta algorithm can be obtained for r = 3 and is determined 
by five independent consistency equations, in the six constants defining the 
algorithm. The equations have multiple solutions, all of which are solutions 
to the equations which determine the classical explicit Runge Kutta algo- 
rithms. Thus, all solutions define classical algorithms, but the solutions are 
not ones traditionally found in the literature. We have used the following 
solution of all five equations: 

= K <^2 = “2/3, C3 = 2/3, 



C21 = -1/24, C3, = 161/24, C32 = -6. 

In the special case of the set of equations (20) and (21), defining the 

elastica on a spliere, we use a hybrid algorithm in which equations (20) are 

integrated using the third order geometrically exact Runge Kutta algorithm 
described above by freezing the coefficients Kg. These coefficients are then 
updated by integrating equation (21) using elliptic functions. 

In equation (20), M is the three dimensional submanifold of deter- 
mined by the equations 

< >= r^, < Xi,X2 >= 0, < Xi.x^ >= 0, 

< X’2,X2 >= 1, < 0, < X3,X3 >= 1. 
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The solutions of (20) for any initial condition .t(0) = xq G M lie completely 
in M. This follows from the fact that the functions 



/, =< a,-i,xi >, /2 =< J-2 >, /a =< >. 

/4 =< X-i, X-i >, /s =< X2, X3>, fe =< X3, X3 > 

satisfy the system of differential equations 

/; = V 2 

f '2 — /4 + ''■5/3 — (1/^^)/] 

/a = /s - ^- 9 h 
n = 2s/ 5 - (2/r^)/2 

n = s/o-(ia-')/3-s/-« 

fe = -2^p/s 

which has the unique solution 

/i = /a = 0, /a = 0, /4 = 1, /s = 0, /e = 1- 

Since this argument does not specif}' the function Kg, it follows that the flows 
of the vector fields F\ — F 3 with frozen coefficients are mappings into M as 
needed in formula (26). 

Furthermore, the vector fields Aj in equation (20) are linear, so that F is 
given by an expression of the form 



F(x) = ^b>(t,x)B,x (27) 

where Bj are matrices and are functions. Freezing the functions to values 
yields a system of linear diflferential equations with constant coefficients: 

i= (i:i}’B,)x. (28) 

^J=I ' 

Thus, the flow Op^ of the vector field Fj (i=l,2,3) is given by 

0p,{t,q) = exp{tCi)q 
and (26) takes the special form 



= exp(c3/tC3) ■ exp(c2/iC2) • exp(c]/tC]) • x^- 
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Since = \/r is constant and = Kg depends only on the matrices C, 
are given by 



C\ — C{ik)t C 2 — C{tk + hc2i), C 3 — C{tk + h{c3i + C 32 )) 

where C{t) ;= C can be considered as a 3 X 3 skew symmetric 

matrix with matrix components that are themselves 3x3 matrices. Therefore 
the standard formula for the exponential of a 3 x 3 skew symmetric matrix 
can be used to compute the flow of each vector field: 

exp(i<;f>S(c)) = 7 + sm{t<j))S{c) + (1 — cos{t(f)))S{c)^ 

where S{c) is the skew symmetric matrix satisfying S{a)b = bx a, and |c| = 1 . 

Regarding the performance of the geometrically stable integration method 
it is clear that this method is computationally more expensive than a classical 
algorithm with the same number of stages and step length. The e.xact cost 
of the new integration scheme can be found in [ 6 ]. However, if we compare 
the performance of the geometrically stable method with the classical fourth 
order Runge-Kutta algorithm on the problem of spherical elastica, it turns 
out that a slightly increased step length for the new method suffices to out- 
perform the classical integration method. Using MATLAB implementations 
of both algorithms we found that the geometrically stable method needs an 
average value of 1001 Mflops for one complete step while the fourth order 
Runge-Kutta only uses an average value of 768 Mflops per step. Therefore, 
one can statistically achieve the same performance by using an increased step 
length of A ~ 1.3038 for the new method. Taking into account that the 
proposed algorithm not only delivers points that lie exactely on the sphere 
but that are also approximately equally spaced along the tracked curve, the 
geometrically stable method seems to be a good choice for the integration of 
spherical elastica. 

A performance comparison of the new integration scheme with a more 
sophisticated classical method, the IMSL Runge Kutta implementation, can 
be found in [7]. 

5 Classification of spherical elastica 

Acting on the suggestion of D. Bernoulli, L. Euler derived differential equa- 
tions for plane elastica and classified the fundamental forms of these curves 
(see [9], [11]). A curvature analysis of the various fundamental cases has been 
given in [3]. 

In this section we classify the forms of spherical elastica based on the 
differential equation ( 21 ) for the geodesic curvature. This equation is of the 
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same form as the equation for tlie curvature of plane elastica and can be 
solved in terms of Jacobi’s elliptic functions in the form: 

dn(Av77i(^ I / ) 

where the positive parameter P of the elliptic function is given by 

2{Kl+2y-a) 

(see [3]). The parameter represents the amplitude of the periodic curva- 
ture function and denotes the value at which K{s,n) = 

To obtain a representation of the curvature in terms of Jacobi’s functions 
with parameter P smaller than 1, one uses the formula 

cn(\/(/c2^ -h2//'2 -a)f2{s - s,n) I (30) 

if a < + 2/7’^. Since the function cn has zeros while dn is positive, the 

above case distinction reflects the main division of elastic curves into those 
where the geodesic curvature changes sign and the other with constant sign 
of their geodesic curvature. 

The change of the forms of spherical elastica while is fixed and a in- 
creases is shown by figures 1-8. The maximum value of the tension parameter 
cr for a real elastic curve on a sphere is according to (29) a = + 2/r^. 

This choice of a corresponds to the dashed circle that is shown in all figures 
for the purpose of orientation. The second curve in figure 1 has a negative 
tension parameter of high absolute value {a = —10000). In comparison to 
figure 2 where a = —30 we observe that decreasing the tension parameter 
has the effect of lowering the amplitude of the curve as it is known from 
the euclidean case (see [3]). The oscillation of the curve in figure 1 is of 
too low amplitude to be visually observed and the curve can not be distin- 
guished from the geodesic determined by the initial conditions. (Note, that 
the geodesic curvature (30) is far from getting flat for a ^ — co but in fact 
approaches a cos function that oscillates with a period that decreases with 
a.) 

A curve with positive tension parameter is shown in figure 3 where a = 
2/r^ = 2. The displayed curve is a special case because the parameter 
l/P = 1/2. Here the geodesic curvature is given by the lemniscate function: 

S(-s) = «m coslemn(K„(.s - 6„,)/2). 

In figures 4 and 5 it is illustrated that with increasing positive a the 
bays of the curves start to overlap until a figure-8-configuration is reached 



(29) 
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Figure 1: Hm = 4, a = —10000 



where all double points of the curve coincide. This happens for a parameter 
cr 7.849. 

While a increases further the curve proceeds through the figure-8-shape 
forming a series of loops with alternating sign of geodesic curvature (see 
figure 6). These loops recede from each other until in the limiting case, when 
cr = 2/r^ + 0.5/cf„, the curve forms a single loop (see figure 7). Here the 
geodesic curvature is given by 

Figure 8 shows that the single loop transforms into a series of loops with 
the same sign of geodesic curvature. With increasing a the loops come closer 
together and finally collapse into a circle when cr = 2/r^ + 
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